Point defect energetics in silicon using the LDA+U method 
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We present the first principles results of point defect energetics in silicon calculated using the 
LDA+U method: a Hubbard type on-site interaction added to the local density approximation 
(LDA). The on-site Coulomb and exchange parameters were tuned to match the experimental band 
gap in Si. The relaxed configuration was obtained using the LDA; LDA+U was used only to 
calculate the energies. Our values of point defect energetics are in very good agreement with both 
recent experimental results and quantum Monte Carlo (QMC) calculations. 

Keywords: Silicon, self-diffusion, point defects, local density approximation, Hubbard, LDA+U 



I. INTRODUCTION 

Modeling techniques empower the investigator with 
the capability to explore into the details of natural 
phenomena with great dexterity, aptly complementing 
their experimental counterparts. Hierarchical multiscale 
modeling (HMM) is one such technique. We have re- 
cently investigated self diffusion in silicon-germanium al- 
loys in considerable detail by developing a database of 
first principles energetics and using the database to per- 
form kinetic Monte Carlo simulations: a typical HMM 
schemei. The exponential growth in computational pro- 
cessing power makes such schemes more viable; the explo- 
sive growth in nanotechnology provides technologically 
relevant platforms that are amenable to such schemes. 

Understandably, the accuracy of the results such as 
those presented in Ref. 0] hinges to a great extent on the 
correctness of the energetics database. In that study, we 
used the popular local density approximation (LDA) to 
develop the energetics database. However, the (approxi- 
mately 1 eV) discrepancy between the theoretical activa- 
tion energy for Si self diffusion computed using the LDA 
(or the generalized gradient approximation (GGA)) and 
the experimental values was an important reason that 
precluded us from making direct comparisons of our re- 
sults with experimental studies of self diffusion in SiGe 
alloys such as those in Refs. 00- 

Our present work is directed essentially at finding a 
means of developing a reliable point defect energetics 
database (like in Ref. Ijj) using an as parameter free 
a technique as possible. (While considerable progress 
has been made on techniques like the quantum Monte 
Carlo (QMC) method, it is beyond the current computa- 
tional capacity to use it to develop a detailed energetics 
database as required for our HMM scheme.) We focus 
on the activation energy for self-diffusion in pure Si. For 
self diffusion in Si, we consider the vacancy mechanism, 
the interstitialcy mechanism, and the concerted exchange 
mechanism. These three have been shown to be the most 



likely diffusion mechanisms in silicon^. For the inter- 
stitialcy mechanism, we consider the hexagonal to split 
(110) mechanism as this has shown to be the most likely 
mechanismi 

This article is organized as follows: In Sec. |H1 we 
provide the details of the LDA and the LDA+U calcu- 
lations. In Sec. IIIII we present the results of our calcu- 
lations. We provide an intuitively appealing explanation 
for the underperformance of the LDA. We then provide 
a brief overview of the LDA+U technique and discuss 
the results obtained using this method for the different 
diffusion mechanisms. We summarize our article in Sec. 

EH 

II. METHOD 

Self-consistent electronic structure calculations were 
performed using the plane-wave pseudopotential code 
VASP 8-9-1011 with the projector augmented-wave (PAW) 
potential o 12 ' 13 at a kinetic energy cutoff of 20 Rydberg. 
A 64 atom supercell was used. Electronic minimization 
was carried out to a tolerance of 2.7 x 10~ 5 eV and struc- 
tures were relaxed until the maximum force on any atom 
was less than 0.015 eV/A. Saddle point configuration 
for the hexagonal to split (110) interstitialcy mechanism 
was determined using the nudged elastic band method 14 . 
Structural relaxation was performed in two stages: ini- 
tially with a 2 3 Monkhorst-Packi^ k-point sampling fol- 
lowed by a 6 3 Monkhorst-Packi^ k-point sampling. En- 
ergy calculations were done using tetrahedron method 
with a 6 x 6 x 6 division of Brillouin zone. The lattice 
constant of systems containing point defects (vacancy or 
interstitial) was determined by fitting the total energy 
versus the supercell volume to Murnaghan's equation of 
stated!. 

The plane- wave psedopotential code VAS P 8 i 9 i 10 i 1:L also 
has the option to perform calculations by the LDA+U 
method. (A brief overview of the method itself will be 
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presented in Sec. IIIII Here we merely give the details 
of the parameters used in our calculations.) We have 
used the rotationally invariant LDA+U scheme according 
to Liechtenstein et ali&. The on-site interactions were 
added for the p-orbitals. The effective on site Coulomb 
interaction parameter (U) was set to eV and the effec- 
tive on site exchange interaction parameter (J) was set to 
4 eV. (We provide justification for this choice in Sec. 11111 ) 
The LDA was used to perform structural relaxation; the 
LDA+U was used only to perform energetics calculation 
on the structure obtained using the LDA. We adopted 
this approach^ because the LDA gives structutal prop- 
erties (eg. lattice constant, see Table that are closer 
to experimental values than the LDA+U. 



III. RESULTS AND DISCUSSION 

Table [I] summarizes the results of both the LDA and 
LDA+U calculations. It also contains experimental val- 
ues and values from quantum Monte Carlo simulations 
where available. 

There is consensus on the values obtained within the 
LDA by different theoretical groups^iii^. It can be seen 
that the LDA gives activation energy that is systemati- 
cally lower than the experimental values. 



TABLE I: Theoretical results and experimental values 
for the following properties [symbol] (units) of silicon: lat- 
tice constant [Lq] (A), vacancy formation volume t [Uy](A 3 ), 
vacancy formation energy[_Ey](eV), vacancy migration 
energy[_E™](eV), interstitial formation volume tt [V/](A 3 ), 
hexagonal interstitial formation energy[_B^ J ](eV), split 
(110) interstitial formation energy [£'g / ](eV), interstitial 
migration energy[i?™](eV), concerted exchange migration 
energy[_EJ? x ](eV). Superscripts indicate reference numbers. 
Quantum Monte Carlo results are from Ref. [2gL 
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* These values are the experimental estimates of the 
activation energy i.e., the sum of formation and migration 
energies. 



A. Local density approximation 

One of the main reasons for the underperformance of 
the LDA is that it tends to overbind. Because of this, the 
dangling bonds that are created due to the defects, tend 
to bind with each other and with other orbitals. Because 
of this overbinding, the energy of the system with the 
defect is lower than it should be in reality. Thus the 
difference between the energy of the pure crystal and the 
one with the defects (which is in fact the formation or the 
migration energy, depending on the situation) is lower 
than it should be in reality. The reason for the LDA's 
overbinding can be traced to the spherical nature of the 
exchange hole. This leads to a larger negative value of 
the exchange energy and hence lowers the energy of the 
system much more than it should be in reality^. 

A related reason for the overbinding nature is the well- 
known underestimation of the band gap by LDA. From 
our LDA based calculations, we find that, for example, 
the Si with a vacancy has a band gap of only 0.88 eV in 
contrast to experimental values of 1.12 eV. (We note that 
the band gap in pure Si from LDA calculations is even 
less: 0.56 eV (Figs. ^HJ- The system with the impurity 
has a higher band gap because of the interaction of the 
defect with the other atoms in the lattice. This can be 
seen from a simple tight binding analysis. We consider 
the band gap of the system with the defect because that 
is the environment that the defect "perceives" .) Because 
of the small band gap, the defect states that are (usually) 
created in the band gap tend to interact more strongly 



with the valence states. This therefore leads to overbind- 
ing too. The band structure and the density of states 
of the system with a vacancy, hexagonal interstitial and 
split (110) interstitial are shown respectively in Figs. |3 
H andEl 

There are a couple of interrelated reasons for the LDA's 
underestimation of the band gap. The discontinuity of 
the one-electron potential for localized states, which is a 
characterisitc of the exact density functional, is absent in 
the LDA as was shown by Perdew et. alM: Because the 
band gap can be expressed as 

E gap = E[N + 1] + E[N - 1] - 2E[N] (1) 

where E[x] is the energy of the system with x electrons, 
the absence of the discontinuity of one-electron potential 
in the LDA causes the LDA to give an incorrect band 
gap. The second reason for the underestimation of band 
gap within LDA is the absence of self interaction correc- 
tion. The third reason has to do with the Kohn-Sham 
approach24 itself. In the Kohn-Sham theory, there is no 
direct relationship between the orbital energies and the 
ionization energies (unlike that in Hartree-Fock theory) 
except for the highest occupied orbital. 

B. Local density approximation with Hubbard 
type correction: LDA+U 

The LDA+U method has been widely applied to metal- 
lic systems containing localized electrons. We refer the 
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reader to the review article on LDA+U by Anisimov et. 
alw^. Here we provide a brief explanation of the tech- 
nique along the lines of "rectifying" the deficiencies of 
the local density approximation that were pointed out 
in Sec. IIII Al . The basic idea of the LDA+U method 
is to redistribute the electron-electron interaction. Bc- 



cause, as we explained in Sec. IIII Al overbinding of the 
localized electrons associated with the dangling bonds is 
a primary reason for the underperformance of the LDA, 
a direct way of rectifying this deficiency would be to add 
an on-site Coulomb repulsion term while subtracting an 
average overall electron-electron interaction. 



X W K 



FIG. 1: Si band structure computed using the LDA (solid 
line) and the LDA+U (dashed line) techniques. The (indirect) 
band gap using LDA (LDA+U) is 0.56 (0.95) eV. 



N N 

Hlda+u = Hlda + Uij / 



where N = J2i=i I IV'i( r )| 2 dr corresponds to the total 
number of electrons in the system. The parameter U is 
adjusted by trial and error so that the resulting band 
gap equals the experimentally observed band gap of the 
system in question. The above technique, which only has 
the Coulomb term, is generalized to inclu de the e xchange 
term as well. We refer the reader to Refs. Ilfil2fil for more 
details. The UN(N — l)/2 introduces the discontinuity 
of the of the one-electron potential. We wish to point out 
that the LDA+U does not correct for the self interaction. 

As mentioned in Sec. [HJ we have used the scheme as 
implemented in VASP for our calculations. Ideally, the 
Coulomb (U) and exchange (J) interaction parameters 
would be chosen so that the band gap for the bulk sys- 
tem would be equal to the experimentally observed band 
gap. However, because our calculations are based on a 
supercell geometry, there is an inevitable but unphysi- 
cal interaction between the defect and the matrix. This 
alters the band gap from that in the bulk. However, 
because the system "sees" only this band gap, we have 
chosen U and J so that the band gap in the system with 
the defect is close to the experimentally observed band 
gap. This, of course, would mean that the parameter 
be tuned for each type of defect. However, for the sake 
of demonstration, we have worked with a single set of 
parameters. 

We find that the defect formation and migration ener- 
gies increase almost linearly with the difference: U — J. 



)| 2 dn J hfc(f2)| 2 dr 2 - U N{N 2 1] (2) 

I 

The best agreement between the above energies and the 
experimental data occurs when the calculated band gap 
corresponds to the experimental one. (This occurs when 
we choose U = and J = AeV as mentioned in Sec. Ull ) 
This supports our argument that the correction of the 
density functional for the bandgap (discontinuity of the 
potential on the number of electrons) should improve the 
results. Table [I] summarizes our energetics calculations. 
Figs. 11151 show band structures and densities of states of 
the various systems considered. 



C. Defect structure, symmetry, states and energy 

1. Vacancy 

There is an inward relaxation of the atoms around the 
vacancy site. The symmetry of the vacancy determines 
the impurity states' symmetry. The singlet s-state is lo- 
cated deep inside the valence band (not shown in Fig. 
|3J) while partly occupied triplet p-states are located in 
the bandgap (FigJSJl. There is a splitting of p- states 
into twofold degenerate lower band and a single upper 
band (almost unoccupied). There are a few major differ- 
ences between the LDA results and the LDA+U results. 
First, there is the expected increase in the bandgap in 
the LDA+U compared to the LDA leading to the separa- 
tion between the vacancy states and valence and conduc- 
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FIG. 2: Band structure and density of states of bulk Si computed using the LDA (left) and the LDA+U (right) methods. This 
figure is for reference to compare with similar figures for systems with defects. 



tion bands. Second, the dispersion of vacancy p-states is 
smaller in the LDA+U calculation. Thus, these states are 
more localized in the LDA+U calculation. The disper- 
sion cannot be removed in the supercell geometry com- 
pletely, but the LDA definitely seems to overestimate the 
itinerancy of electrons in the vacancy states. There is no 
noticeable Jahn- Teller distortion of the atoms surround- 
ing the vacancy in our 64-atom supercell geometry. This 
is consistent with the observation by Zywietz et. alwL 
who suggest that a 128-atom supercell is required to ob- 
serve the distortion. The formation and the migration 
energy for the vacancy defect from LDA+U calculations 
are 4.59 eV and 0.42 eV respectively. This is in good 
agreement with experimental results^S^ 



2. Hexagonal interstitial 

The self interstitial sites provide four additional elec- 
trons to the system and create similar "dangling" bonds. 
Like the vacancy, they create an almost doubly degen- 
erate band near the bottom of the band gap and a sin- 
gle band near the top of the band gap (Fig. 0}. The 
occupied defect states are very close to the bottom of 
the band gap, in agreement with the observation made 
by Needsi about the defect states being quite shallow. 
We find that the atoms surrounding the interstitial move 
outwards in contrast to the inward movement reported 
by Needst This is probably because, while we have op- 
timized the lattice constant of the system with the de- 
fect, Needs has maintained the lattice constant at the 
experimental value. We suspect that Needs' approach 
might cause an unphysical restriction on the relaxation. 
In addition to opening the band gap, the introduction 
of the on-site repulsion quite dramatically reduces the 
mixing of the defect states with the valence states. The 



formation energy for the hexagonal interstitial using the 
LDA+U method is 4.61 eV which is in good agreement 
with experimental- and quantum Monte CarloSi results. 



3. Split (110) interstitial 

The split (110) interstitial creates defect states in the 
band gap similar to the hexagonal interstitial excpet that 
there is a greater splitting of the two lower bands (Fig. 
|SJ). We attribute this to the lower symmetry of the split 
(110) interstitial compared to the hexagonal interstitial. 
Another difference is the location of the higher defect 
state. It is not as close to the top as in the case of hexag- 
onal interstitial. Wc do not have a simple explanation 
for this. The introduction of the on-site repulsion has 
a similar effect with respect to reducing the mixing of 
the defect states with the bulk states. The defect forma- 
tion energy using the LDA+U method is 4.70 eV which 
is in good agreement with the experimental and quan- 
tum Monte Carlo^S results. The migration energy from 
the hexagonal to the split (110) configurations of the in- 
terstitial was computed using the nudged elastic band 
menthodii We get a value of 0.44 eV using the LDA+U 
method. We are not aware of any experimental result for 
this specific value. 



D. Self diffusion in Si 

The contribution of the mechanism to the diffusion 
process is determined by the activation energy and the 
prefactor. Because all the three mechanisms considered 
in this article (viz., vacancy, interstitial, and concerted 
exchange) involve similar number of atoms, one can, to 
a first approximation, assume that all of them have sim- 
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FIG. 3: Band structure and density of states of a system containing a vacancy computed using the LDA (left) and the LDA+U 
(right) methods. The defect states are shown as dashed lines. The (indirect) band gap using the LDA (LDA+U) is 0.88 (1.21) 
eV and the vacancy dispersion using the LDA (LDA+U) is 0.77 (0.71) eV. The bottom vacancy state is doubly degenerate. 
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FIG. 4: Band structure and density of states of a system containing a hexagonal interstitial computed using the LDA (left) and 
the LDA+U (right) methods. The defect states are shown as dashed lines. The (indirect) band gap using the LDA (LDA+U) 
is 0.70 (1.01) eV. The amount of mixing of the defect states with the valence states drops in the LDA+U method, signifying 
lesser binding. 



ilar entropic effects and hence similar prefactors. Thus, 
based on our calculations of the activation energies using 
the LDA+U method, the vacancy mechanism [which has 
an activation energy of 5.01 eV (sum of formation (4.59 
eV) and migration (0.42 eV) energies)] and the hexag- 
onal to split (110) interstitialcy mechanism [which has 
an activation energy of 5.14 eV (sum of formation (4.7 
eV) and migration (0.44 eV) energies)] are quite close in 
activation energy (within the error of the method), and 
should contribute equally to self-diffusion. The concerted 
exchange mechanism (which has an activation energy of 
5.82 eV) is less significant in this respect. Experimental 
results give a similar indication 6,29 . 



IV. SUMMARY 



The present work identified the cause for the poor 
description of point defect energetics by the LDA. It 
corrected for the deficiencies of the LDA by using the 
LDA+U method. This gave much better agreement of 
the calculated activation energies with experimental ob- 
servations. This method can therefore be used for better 
description of diffusion in similar semiconductor materi- 
als. 
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FIG. 5: Band structure and density of states of a system containing a split (110) interstitial computed using the LDA (left) and 
the LDA+U (right) methods. The defect states are shown as dashed lines. The (indirect) band gap using the LDA (LDA+U) 
is 0.99 (1.12) eV. The amount of mixing of the defect states with the valence states drops in the LDA+U method signifying 
lesser binding. 
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